Linear Regression¶
Setup¶
In [277]:
import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas_profiling
from sklearn.datasets import load_boston
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import ElasticNet
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from sklearn.pipeline import FeatureUnion
from sklearn.pipeline import Pipeline
from sklearn.base import BaseEstimator
from sklearn.base import TransformerMixin
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import OneHotEncoder
from sklearn.model_selection import train_test_split
#conda install -c districtdatalabs yellowbrick --yes
from yellowbrick.regressor import ResidualsPlot
from yellowbrick.regressor import PredictionError
import statsmodels.api as sm
import statsmodels.stats.api as sms
import statsmodels.formula.api as smf
import statsmodels.graphics.regressionplots as regplot
from statsmodels.stats.outliers_influence import variance_inflation_factor
from patsy import dmatrices
Data¶
In [278]:
boston = load_boston()
print(boston['DESCR'])
In [279]:
bostonDF = pd.DataFrame(boston['data'], columns=boston['feature_names'])
bostonDF['TARGET'] = boston['target']
bostonDF.head()
Out[279]:
Categorical Features¶
- We'll fabricate age into a catagorical feature
In [280]:
def ageBin(age):
if age < 30: return 'NEW'
if age < 60: return 'USED'
if age < 80: return 'OLD'
if age < 95: return 'CREAKY'
return 'FOSSIL'
bostonDF['AGE_BIN'] = bostonDF['AGE'].apply(ageBin)
bostonDF['AGE_BIN'].hist()
Out[280]:
EDA¶
In [281]:
pandas_profiling.ProfileReport(bostonDF)
Out[281]:
Preprocessing¶
Flow Helpers¶
In [282]:
class ColumnSelector(BaseEstimator, TransformerMixin):
def __init__(self, columns):
self.columns = columns
def fit(self, X, y=None):
return self
def transform(self, X):
assert isinstance(X, pd.DataFrame)
try:
return X[self.columns]
except KeyError:
cols_error = list(set(self.columns) - set(X.columns))
raise KeyError("The DataFrame does not include the columns: %s" % cols_error)
class TypeSelector(BaseEstimator, TransformerMixin):
def __init__(self, dtype):
self.dtype = dtype
def fit(self, X, y=None):
return self
def transform(self, X):
assert isinstance(X, pd.DataFrame)
return X.select_dtypes(include=[self.dtype])
Age¶
In [283]:
bostonDF['AGE_BIN'] = bostonDF['AGE_BIN'].astype('category')
bostonDF = bostonDF.drop('AGE', axis=1)
Build Pipeline¶
In [284]:
process = make_pipeline(
FeatureUnion(transformer_list=[
("numeric_features", make_pipeline(
TypeSelector(np.number),
StandardScaler()
)),
("categorical_features", make_pipeline(
TypeSelector("category"),
OneHotEncoder(categories='auto', drop='first')
))
])
)
Loss Metrics¶
Explained Variance¶
- Talking point - What is R-Squared?
- ESS / TSS
MSE¶
- Measures the average squared error magnitude
- Not in units of the prediction
- Hard to explain / conceptualize
RMSE¶
- In the same units of the prediction
- Measures the average magnitude of the error
- Larger errors result in larger RMSE (from squaring the term)
MAPE & Accuracy¶
Model¶
Tracking¶
In [285]:
resultsDF = pd.DataFrame(columns=['MAE', 'MSE', 'RMSE'])
compareDF = pd.DataFrame()
Fit¶
In [286]:
df = bostonDF.copy()
y = df['TARGET']
X = df.drop('TARGET', axis=1)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1337)
clf = LinearRegression()
pipe = Pipeline(steps=[('process', process), ('regression', clf)])
pipe.fit(X_train, y_train)
pred = pipe.predict(X_test)
mae = mean_absolute_error(y_test, pred)
mse = mean_squared_error(y_test, pred)
rmse = math.sqrt(mse)
print('MAE = {:.2f}'.format(mae))
print('MSE = {:.2f}'.format(mse))
print('RMSE = {:.2f}'.format(rmse))
resultsDF.loc['OLS'] = [mae, mse, rmse]
compareDF = pd.concat([compareDF, pd.DataFrame(clf.coef_)], axis=1)
In [324]:
errors = abs(y_test - pred)
mape = 100 * (errors / y_test)
accuracy = 100 - np.mean(mape)
print('MAPE : {:.2f}%'.format(np.mean(mape)))
print('Accuracy : {:.2f}%'.format(accuracy))
Prediction Error¶
In [287]:
plt.figure(figsize=(15,10))
visualizer = PredictionError(pipe)
visualizer.fit(X_train, y_train) # Fit the training data to the visualizer
visualizer.score(X_test, y_test) # Evaluate the model on the test data
g = visualizer.poof() # Draw/show/poof the data
Residuals Plot¶
- We should look for a cloud like shape (no other pattern)
In [288]:
plt.figure(figsize=(15,10))
visualizer = ResidualsPlot(pipe)
visualizer.fit(X_train, y_train) # Fit the training data to the model
visualizer.score(X_test, y_test) # Evaluate the model on the test data
visualizer.poof()
Interpretation¶
StatsModel¶
In [289]:
tDF = pd.concat([X_train, y_train], axis=1)
f = 'TARGET ~ CRIM + ZN + INDUS + CHAS + NOX + RM + DIS + RAD + TAX + PTRATIO + B + LSTAT + AGE_BIN'
lm = smf.ols(formula=f, data=tDF).fit()
Summary¶
In [290]:
print(lm.summary())
Partial Regression Plots¶
Creates a partial regression subplot for each explanatory variable. Each plot shows the relationship between the response and the given explanatory variable after removing the effect of all other explanatory variables. These plots are helpful in describing the effect of one variable's effect.
In [291]:
fig = plt.figure(figsize=(12,30))
fig = sm.graphics.plot_partregress_grid(lm, fig=fig)
Assumptions¶
L (inear)¶
- Data should be visibly linear - hard to do in high dimensional spaces
- You can attempt to fit polynomial independent variables or apply other transformation to capture / fit the non-linerity
I (ndependent)¶
- Each data point should be independent.
- We want to avoid colinearity
- When we have high correlation on two features, there is little information from both in the model
- Consider dropping out one of the correlated features
- Example violations of this assumption
- Include time series (which can introduce autocorrelation)
- Height regression that includes height of left leg and height of right leg in model
Variable Inflation Factor¶
Variance Inflation Factor (VIF) measures the colinearity of predictor variables. VIF is the ratio of the variance of all model's coefficients divided by the variance of a single coefficient if it were fit in isolation.
In [292]:
variables = lm.model.exog
vif = [variance_inflation_factor(variables, i) for i in range(variables.shape[1])]
vifDF = pd.DataFrame()
vifDF['Coefficient'] = lm.params
vifDF['VIF'] = vif
vifDF['Flag'] = np.where(vifDF.VIF > 5, True, False)
vifDF = vifDF.drop('Intercept', axis=0)
vifDF
Out[292]:
Correlation RAD & TAX¶
- These two features are highly correlated and score high on the VIF.
- RAD = index of accessibility to radial highways
- TAX = full-value property-tax rate per $10,000
In [293]:
print(np.corrcoef(df.TAX, df.RAD))
N (ormally Disributed Residuals)¶
- Residuals should be cloud shaped with no discernable pattern
- Check highly influential points (possible outliers)
Residual Plot¶
In [294]:
preds = lm.predict(X_train)
plt.scatter(preds, -lm.resid)
plt.title('Residual Plot (Residuals vs Pred)')
plt.ylabel('Residual')
Out[294]:
Jarque Bera Test for Normality¶
In [295]:
cols = ['Jarque-Bera', 'Chi^2 two-tail prob.', 'Skew', 'Kurtosis']
test = sms.jarque_bera(lm.resid)
pd.DataFrame(zip(cols, test)).T
Out[295]:
Leverage & Influence Plot¶
Leverage
- The ability to exert influence over the model
- Measures of how far an observation on the predictor variable is from from the mean of the predictor.
- The higher the leverage value, the more potential it has to impact the fitted model.
Influence
- Degree to which the data point affects the model.
- Measures the effect of a point on the regression model's predictions.
- A high influence point usually has leverage.
Watch for:
- High leverage and large residual - No points
- High Leverage and small residual - Point (380)
- Low Leverage and large residual - Points (372, 368, 364)
In [296]:
fig, ax = plt.subplots(figsize=(12,8))
fig = sm.graphics.influence_plot(lm, ax=ax, criterion="cooks")
High Leverage Points¶
In [297]:
points = [380, 418]
X_train[X_train.index.isin(points)]
Out[297]:
In [298]:
pred = lm.predict(X_train)
plt.title('High Leverage Points (Predicted vs Actual)')
plt.scatter(y_train, pred)
plt.ylabel('Actual')
for pt in points:
plt.scatter(pred.loc[pt], y_train.loc[pt], color='red')
In [299]:
preds = lm.predict(X_train)
plt.scatter(preds, -lm.resid)
plt.ylabel('Residual')
plt.title('Predicted vs Residual')
for pt in points:
plt.scatter(pred.loc[pt], -lm.resid[pt], color='red')
Large Residual Points¶
In [300]:
points = [364, 372, 368]
X_train[X_train.index.isin(points)]
Out[300]:
In [301]:
pred = lm.predict(X_train)
plt.title('Large Residual Points (Predicted vs Actual)')
plt.scatter(y_train, pred)
plt.ylabel('Actual')
for pt in points:
plt.scatter(pred.loc[pt], y_train.loc[pt], color='red')
In [302]:
preds = lm.predict(X_train)
plt.scatter(preds, -lm.resid)
plt.ylabel('Residual')
plt.title('Predicted vs Residual')
for pt in points:
plt.scatter(pred.loc[pt], -lm.resid[pt], color='red')
E (rror constant)¶
- Homoskedasticity vs Heteroskedasticity
- Example - As income goes up our predicted error increases
- Example - As we play more minutes of basketball our predicted scoring error increases
- Typically spot with visual plot of residual vs predicted value (look for cone shape)
- If find heteroskedasticity we typically:
- Adjust the model - add / drop featurews
- Box Cox (log / sqrt / power) transforms on features and the target
Visual Inspection¶
- There is some structure to the residual plot.
- Indicative of heteroksedaticity
In [303]:
preds = lm.predict(X_train)
plt.scatter(preds, -lm.resid)
plt.ylabel('Residual')
plt.title('Predicted vs Residual')
plt.show()
Breush Pagan Test¶
Produces a p-value that when less that a significance level of 0.05, therefore we can reject the null hypothesis that the variance of the residuals is constant and infer that heteroscedasticity is present.
In [304]:
cols = ['Lagrange multiplier statistic', 'p-value', 'f-value', 'f p-value']
test = sms.het_breuschpagan(lm.resid, lm.model.exog)
pd.DataFrame(zip(cols, test)).T
Out[304]:
Regularization¶
Adding features to Linear Regression will always have the effect of improving R-Squared. To mitigate the chance of overfitting (fitting noise in the data) we want to penalize higher larger coefficients using L1 or L2 Norms. The choice of L1 or L2 has an important effect.
Norms¶
- Norm is the length of all vectors in a space
- L1 Norm
- Manhattan Distance
- Sum of absolute values of coefficients
- L2 Norm
- Euclidean Distance
- Sum of squared coefficients
- L1 and / or L2 are added to the cost function
- Large coefficients result in higher penalty
Lasso - L1 Regularization¶
Model¶
In [305]:
df = bostonDF.copy()
y = df['TARGET']
X = df.drop('TARGET', axis=1)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1337)
clf = Lasso()
pipe = Pipeline(steps=[('process', process), ('regression', clf)])
pipe.fit(X_train, y_train)
pred = pipe.predict(X_test)
mae = mean_absolute_error(y_test, pred)
mse = mean_squared_error(y_test, pred)
rmse = math.sqrt(mse)
print('MAE = {:.2f}'.format(mae))
print('MSE = {:.2f}'.format(mse))
print('RMSE = {:.2f}'.format(rmse))
resultsDF.loc['Lasso'] = [mae, mse, rmse]
compareDF = pd.concat([compareDF, pd.DataFrame(clf.coef_)], axis=1)
Alpha Hyperparameter¶
In [306]:
hyperDF = pd.DataFrame()
alpha = [0.0, 0.05, 0.25, 0.5, 0.75, 1.0]
cols = []
for a in alpha:
col = 'Alpha-{}'.format(a)
clf = Lasso(alpha=a)
pipe = Pipeline(steps=[('process', process), ('regression', clf)])
pipe.fit(X_train, y_train)
hyperDF[col] = clf.coef_
hyperDF = hyperDF.sort_values('Alpha-0.0')
In [307]:
f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12,8), sharey=True)
ax1 = hyperDF['Alpha-0.0'].plot.bar(ax=ax1)
ax1.set_title('Lasso - Alpha 0.0')
ax2 = hyperDF['Alpha-0.25'].plot.bar(ax=ax2)
ax2.set_title('Lasso - Alpha 0.25')
ax3 = hyperDF['Alpha-1.0'].plot.bar(ax=ax3)
ax3.set_title('Lasso - Alpha 1.0')
Out[307]:
Ridge - L2 Regularization¶
Model¶
In [308]:
df = bostonDF.copy()
y = df['TARGET']
X = df.drop('TARGET', axis=1)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1337)
clf = Ridge()
pipe = Pipeline(steps=[('process', process), ('regression', clf)])
pipe.fit(X_train, y_train)
pred = pipe.predict(X_test)
mae = mean_absolute_error(y_test, pred)
mse = mean_squared_error(y_test, pred)
rmse = math.sqrt(mse)
print('MAE = {:.2f}'.format(mae))
print('MSE = {:.2f}'.format(mse))
print('RMSE = {:.2f}'.format(rmse))
resultsDF.loc['Ridge'] = [mae, mse, rmse]
compareDF = pd.concat([compareDF, pd.DataFrame(clf.coef_)], axis=1)
Alpha Hyperparameter¶
In [320]:
hyperDF = pd.DataFrame()
alpha = [0.0, 0.05, 0.25, 0.5, 0.75, 1.0, 500.0]
cols = []
for a in alpha:
col = 'Ridge-{}'.format(a)
clf = Ridge(alpha=a)
pipe = Pipeline(steps=[('process', process), ('regression', clf)])
pipe.fit(X_train, y_train)
hyperDF[col] = clf.coef_
hyperDF = hyperDF.sort_values('Ridge-0.0')
In [321]:
f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12,8), sharey=True)
ax1 = hyperDF['Ridge-0.0'].plot.bar(ax=ax1)
ax1.set_title('Ridge - Alpha 0.0')
ax2 = hyperDF['Ridge-1.0'].plot.bar(ax=ax2)
ax2.set_title('Ridge - Alpha 1.0')
ax3 = hyperDF['Ridge-500.0'].plot.bar(ax=ax3)
ax3.set_title('Ridge - Alpha 500.0')
Out[321]:
Ridge vs Lasso¶
- Lasso moves several exogenous variables to 0.0
In [322]:
fig, ax = plt.subplots(figsize=(10,8))
compareDF.columns = ['OLS', 'Lasso', 'Ridge']
compareDF = compareDF.sort_values('OLS')
plt.title('Comparison of Model Components')
plt.ylabel('Coefficient Value')
compareDF[['OLS', 'Ridge', 'Lasso']].plot.bar(ax=ax)
Out[322]:
ElasticNet - Lasso + Ridge¶
In [312]:
df = bostonDF.copy()
y = df['TARGET']
X = df.drop('TARGET', axis=1)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1337)
for l1 in np.arange(0.1, 1, 0.1):
clf = ElasticNet(l1_ratio=l1)
pipe = Pipeline(steps=[('process', process), ('regression', clf)])
pipe.fit(X_train, y_train)
pred = pipe.predict(X_test)
mae = mean_absolute_error(y_test, pred)
mse = mean_squared_error(y_test, pred)
rmse = math.sqrt(mse)
classifier = 'ElasticNet - L1={:.1f} L2={:.1f}'.format(l1, 1.0 - l1)
print ('{} MAE {:.3f} MSE {:.3f} RMSE {:.3f}'.format(classifier, mae, mse, rmse))
resultsDF.loc[classifier] = [mae, mse, rmse]
Comparison¶
- ElasticNet L1 = 1.0 == Lasso
- ElasticNet L2 = 1.0 == Ridge (Warning unsable l1_ratio < 0.01)
- Lasso or Ridge with Alpha = 0.0 == OLS
Top By RMSE¶
In [313]:
resultsDF.sort_values('RMSE')
Out[313]:
Top By MAE¶
In [314]:
resultsDF.sort_values('MAE')
Out[314]:
Notes¶
- Consider regularization as a default for predicting Linear Models
- Ridge promotes coefficients close to zero
- Lasso promotes sparsity and can be used for feature selection
- In this example, Lasso kept 3 - 4 non-zero features + Intercept
- Ridge kept all features
- Lasso / Ridge can change coefficients quite dramatically (positive to negative) with different levels of alpha
- Remember Assumptions
- LINE
- Inspect Residuals & Prediction Plots
- Look for High Leverage Points
- Produce multiple metrics
- Remember differences between MAE & RMSE
- R-Squared is intuitive but limited
- Stick to SKlearn if your objective is highest accuracy.
- Easy to swap in / out estimators
- Pipeline support is helpful
- Consider Statsmodel / R if you're modeling Linear Regression and want deep explanatory capabaility
- Statsmodel has support for Lasso / Risge with .fit_regularized (0 <= L1_wt <= 1)